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FOREWORD 


The  Defense  Communications  Engineering  Center  (DCEC)  Technical  Notes 
(TN's)  are  published  to  inform  interested  members  of  the  defense  community 
regarding  technical  activities  of  the  Center,  completed  and  in  progress.  They 
are  intended  to  stimulate  thinking  and  encourage  information  exchange;  but 
they  do  not  represent  an  approved  position  or  policy  of  DCEC,  and  should  not 
be  used  as  authoritative  guidance  for  related  planning  and/or  further  action. 

Comments  or  technical  inquiries  concerning  this  document  are  welcome,  and 
should  be  directed  to: 

Director 

Defense  Comnunlcations  Engineering  Center 
1860  Wiehle  Avenue 
Reston,  Virginia  22090 


EXECUTIVE  SUMMARY 


Survivability  of  uninterrupted  slip-free  communications  in  a  synchronous 
switched  digital  network  depends  upon  the  survivability  of  the  timing 
function.  The  survivability  of  the  timing  function  can  be  greatly  enhanced  by 
providing  a  backup  free-running  clock  mode  of  timing  for  each  node.  The 
usefulness  of  a  clock  in  the  free-running  mode  is  dependent  on  its  stability 
and  accuracy.  This  becomes  more  Important  at  higher  data  rates  such  as  those 
for  the  future  DCS.  For  stability  and  accuracy  the  DSCS  uses  cesium  clocks 
which  are  controlled  to  keep  them  within  a  few  microseconds  of  the  Naval 
Observatory  master  clock.  Cesium  clocks  could  also  be  used  in  the  rest  of  the 
DCS,  but  they  are  very  expensive.  Other  means  exist  for  providing  more 
accuracy  and  stability  during  a  controlled  mode  of  operation  and  for  short 
periods  of  time  during  free-running  modes.  These  other  means  are  also  more 
survivable  than  the  means  presently  used  to  coordinate  the  DSCS  clocks  with 
the  Naval  Observatory  master  clock.  There  is  an  excellent  possibility  that 
high-quality  quartz  crystal  or  rubidium  clocks  could  be  used  as  a  much  lower 
cost  alternative  for  expensive  cesium  clocks  in  many  DCS  applications  provided 
their  errors  in  a  free-running  period  could  be  predicted  and  removed  with 
sufficient  accuracy  for  a  long  enough  period  of  time.  A  study  by  the  Naval 
Observatory  under  DCA  sponsorship  shows  that  prediction  techniques  can  greatly 
improve  the  accuracy  and  stability  of  clocks  during  a  free-running  period 
following  a  period  of  controlled  operation.  During  the  period  of  controlled 
operation,  referred  to  here  as  a  calibration  period,  errors  in  the  clock  are 
accurately  measured.  These  measured  errors  are  used  to  establish  initial - 
conditions  for  a  mathematical  model  which  is  used  to  predict  clock  errors 
during  the  free-running  period.  The  results  of  the  study  show  that  there  is 
an  extremely  high  probability  of  successful  application  of  the  technique  in 
the  DCS.  This  offers  a  high  potential  capability  for  both  improving 
survivability  and  reducing  costs  compared  to  other  alternatives.  However, 
more  study  of  clock  predictability  using  a  larger  number  of  newer  clocks  is 
needed.  Optimization  of  the  mathematical  models  should  be  attempted,  and 
evaluation  of  possible  aliasing  due  to  infrequent  measurements  (one  per  hour) 
should  be  evaluated.  In  addition  to  these  evaluations  of  the  predictability 
of  clock  performance,  the  study  needs  to  be  extended  to  include  evaluation  of 
the  practical  application  of  the  technology  in  the  DCS. 

This  technical  note  provides  a  background  of  the  need  for  accurate,  stable 
timing  in  a  military  switched  digital  communications  system  such  as  the  DCS. 

It  explains  the  application  of  ARIMA  models  to  clock  prediction,  and  it 
discusses  many  things  that  need  to  be  done  to  answer  questions  about  the 
practical  application  of  these  techniques  to  the  DCS  timing  system. 
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I.  INTRODUCTION 


This  technical  note  presents  a  proposed  program  of  test  and  analysis  with 
a  goal  of  using  prediction  techniques  to  enhance  the  performance  of 
quartz-crystal  clocks,  rubidium  gas  cell  clocks,  and  cesium  beam  clocks  for 
application  in  the  Defense  Communications  System  (DCS).  Although  this  is 
basically  a  research  project,  previous  work  sponsored  by  DCA  indicates  a  very 
high  probability  of  providing  results  which  will  make  it  possible  to  enhance 
the  survivability  of  the  timing  function  or  reduce  costs  or  both.  The 
prediction  models  to  be  employed  are  members  of  the  general  class  of 
Autoregressive  Integrated  Moving  Average  (ARIMA)  models.  Because  most  readers 
are  probably  not  familiar  with  ARIMA  models,  a  discussion  of  these  models  as 
they  can  be  applied  to  clock  prediction  is  also  included. 

Studies  of  timing  (synchronization)  for  the  digital  Defense  Communications 
System  (DCS)  have  shown  the  importance  of  having  a  backup  free-running  clock 
mode  of  -pperation  to  enhante  communications  system  survivability  by  providing 
timing  when  the  normal  timing  capability  is  not  available  for  any  reason.  The 
usefulness  of  a  clock  in  the  free-running  mode  is  dependent  on  its  stability 
and  accuracy.  In  addition  to  reducing  the  need  for  resynchronization  or  the 
occurrence  of  other  communication  outages,  high  stability  and  accuracy  can  aid 
in  system  monitoring  and  can  be  a  most  important  factor  in  providing  rapid 
recovery  of  the  communications  system  following  an  outage. 

The  Defense  Satellite  Communications  System  (DSCS)  presently  uses 
expensive  cesium  clocks  to  aid  the  rapid  acquisition  of  its  spread  spectrum 
signals.  These  clocks  are  coordinated  with  the  Naval  Observatory  master 
clock.  The  techniques  used  for  this  coordination  do  not  have  a  high  degree  of 
survivability,  but  when  adequately  modified  to  provide  the  needed 
survivability,  they  would  also  be  effective  in  other  parts  of  the  DCS. 

However,  cesium  clocks  are  very  expensive  (approximately  $40,000  depending  on 
options),  and  a  less  expensive  alternative  is  needed.  A  preliminary  study  by 
the  Naval  Observatory  under  DCA  sponsorship  shows  a  capability  to  predict 
errors  of  high  quality  quartz-crystal  oscillators  ($8,000  or  less  depending  on 
type  selection  and  options)  or  rubidium  clocks  ($21,000  or  less  depending  on 
type  and  options)  with  considerable  accuracy  during  a  free-running  period 
following  a  calibration  period.  Using  these  predictions,  there  is  a  very  high 
probability  that  errors  can  be  removed  with  sufficient  accuracy  for  a  long 
enough  period  of  time  to  permit  the  use  of  these  clocks  as  an  alternative  to 
much  more  expensive  cesium  clocks  for  many  communications  applications.  An 
important  initial  applies*'  ...  might  use  them  as  alternate  or  backup  clocks  for 
the  principal  cesium  clocks  at  some  DSCS  stations.  Later,  much  more  extensive 
applications  would  be  expected  in  other  parts  of  the  DCS. 

An  objective  of  this  technical  note  is  to  provide  sufficient  background  on 
the  timing  needs  of  the  future  digital  DCS  and  methods  of  satisfying  them  so 
that  the  reader  will  understand:  (1)  what  is  needed  from  this  program  to 
enhance  the  performance  of  clocks  for  application  in  the  DCS,  (2)  why  it  is 
important,  and  (3)  how  it  fits  into  a  complete  timing  system. 

This  document  describes  ARIMA  prediction  models  and  how  they  can  be 
applied  to  the  DCS  so  that  the  reader  can  understand  what  the  models 
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accomplish,  how  they  accomplish  it,  and  the  high  probability  that  their 
application  in  the  DCS  could  result  in  very  significant  cost  savings  and 
enhancement  of  survivability. 

An  objective  of  the  study  program  described  in  this  technical  note  is  to 
further  refine  clock  prediction  models  for  three  types  of  clocks:  (1)  high 
quality  quartz-crystal  clocks,  (2)  rubidium  gas  cell  clocks,  and  (3)  cesium 
beam  clocks.  A  further  objective  is  to  evaluate  the  resulting  models  using 
clock  error  data  obtained  by  comparing  each  of  a  number  of  different  clocks, 
from  each  of  the  three  generic  types,  with  the  Naval  Observatory  master 
clock.  The  evaluations  are  to  be  made  over  a  large  sample  of  different 
periods  of  time  in  order  to  obtain  statistically  meaningful  results. 

Another  aim  of  the  study  program  is  the  development  of  simple  recursive 
microprocessor  algorithms  for  both  calibration  and  prediction  (hopefully,  ones 
which  can  also  permit  the  optimization  of  parameters)  which  will  permit 
practical  application  of  the  technology  to  a  digital  DCS. 

A  further  objective  of  the  study  program,  which  hopefully,  will  also 
include  a  capability  for  optimizing  the  parameters  in  the  prediction  models, 
is  to  include  in  the  programming  of  the  microprocessor  the  capability  to 
filter  timing  information  and  to  control  the  timing  system.  This  includes: 

•  Elimination  of  obviously  erroneous  perturbed  timing  measurements. 

•  Filtering  the  timing  information  during  the  controlled  mode  of 
operation  (most  of  the  time). 

•  Control  of  a  micro-phase-stepper  that  controls  the  output  of  the 
oscillator  or  clock  during  the  controlled  mode  of  operation. 

e  Calibration  of  the  clock  prediction  function  during  the  controlled  mode 
of  operation. 

•  Determination  of  when  it  is  necessary  to  change  from  the  controlled 
mode  of  operation  to  the  free-running  mode. 

•  Making  a  smooth  transition  from  the  controlled  mode  to  the  free-running 
mode. 

•  Performance  of  clock  correction  while  in  the  free-running  mode  (control 
of  the  micro-phase-stepper  to  remove  predicted  errors). 

•  Determination  of  when  to  return  to  the  controlled  mode. 

e  Control  of  a  smooth  transition  back  to  a  controlled  mode  of  operation. 


•  Accomplishment  of  any  other  timing  system  control  functions  discussed 
in  the  references. 


II  BACKGROUND 


As  the  DCS  becomes  more  digital — digital  transmission,  digital  switching, 
digital  control — the  importance  of  system  timing  increases.  The  timing 
relationship  determines  to  whom  the  information  belongs  and  what  it  means; 
i.e.,  particular  time  slots  are  assigned  for  particular  purposes.  It  is  very 
important  in  a  synchronous  digital  network  that  any  bit  originating  anywhere 
in  the  network  be  available  at  the  instant  when  it  is  needed  at  any  node 
through  which  it  passes.  The  loss  of  satisfactory  timing  can  be  catastrophic, 
causing  all  received  information  to  be  meaningless.  Note  that  this  is  a  phase 
control  problem  and  not  just  a  frequency  control  problem,  i.e.;  syntonization 
is  not  adequate  and  synchronization  is  required. 

Studies  by  DCA  and  its  contractors  over  a  number  of  years,  using  both 
simulation  and  analysis,  have  recommended  that  all  major  nodes  of  the  DCS  be 
referenced  to  Coordinated  Universal  Time  (UTC)  when  it  can  be  made  available. 
(References  [1-11]).  When  UTC  is  not  available  to  the  network,  it  is 
recommended  that  all  major  nodes  be  referenced  to  the  best  (highest  ranking) 
clock  in  the  network,  and  to  provide  a  similar  capability  for  any  portion  of 
the  network  that  becomes  isolated  from  the  rest  of  the  network.  These  studies 
recommended  that  a  time  reference  for  all  major  nodes  be  distributed  through 
the  network  by  coordinating  the  transmission  of  synchronization  codes  from 
every  major  node  with  the  network  reference  clock;  that  minor  nodes  be  slaved 
to  major  nodes;  and  that  a  free-running  clock  mode  of  operation  be  provided 
for  each  node,  to  be  used  when  other  modes  of  operation  are  not  available. 
Although  reliability  was  a  concern  in  those  studies,  the  major  concern  has 
more  recently  been  survivability  during  and  following  an  attack  upon  the 
communications  system  prior  to  or  during  full  scale  war.  This  has  led  to 
development  of  a  number  of  attributes  for  timing  in  a  digital  DCS.  A 
discussion  of  these  attributes  [12]  indicates  their  importance  to  a  survivable 
slip-free  digital  DCS.  They  can  be  provided  by  methods  described  in  [13]. 

When  these  attributes  are  provided,  an  accurate  timing  reference  is  available 
at  each  major  node.  This  timing  reference  can  be  used  to  accurately  determine 
the  errors  in  the  local  clock. 

It  was  early  recognized  that  the  stability  and  accuracy  of  a  timing  system 
could  contribute  greatly  to  the  ability  of  a  node  to  reenter  the  network 
quickly  after  a  communications  outage  during  which  its  clock  must  free-run. 
Because  of  this,  stable  and  expensive  cesium  clocks  which  are  maintained  by 
the  Naval  Observatory  within  acceptable  tolerances  of  UTC  (USNO)  have  been  in 
operation  in  the  Defense  Satellite  Communications  System  (DSCS)  for  a  number 
of  years.  By  accurately  distributing  time  to  major  nodes  throughout  the 
remainder  of  the  DCS,  a  very  stable  network  timing  system  will  be  provided. 
This  will  also  provide  the  advantages  of  accurate  clocks  at  all  nodes  which 
can  be  used  for  other  purposes.  With  the  ability  to  accurately  measure  errors 
in  local  clocks  at  major  nodes  of  the  DCS,  it  seemed  there  would  be  an 
excellent  opDortunity  to  use  this  information  to  predict  clock  errors  which 
would  occur  during  a  free-running  period  immediately  following  a  period  of 
measurement  and  calibration.  These  predicted  errors  could  then  be  removed 
during  the  free-running  period,  greatly  improving  the  clock  performance.  This 


might  make  it  possible  to  replace  some  expensive  cesium  clocks  with  much  less 
expensive,  high  quality  quartz-crystal  clocks.  For  those  applications  where 
errors  in  quartz-crystal  clocks  cannot  be  predicted  with  sufficient  accuracy 
over  a  sufficiently  long  period  of  time,  it  is  quite  likely  that  errors  in 
much  more  stable  rubidium  clocks  can  be  adequately  predicted.  The  cost  of 
rubidium  clocks,  though  higher  than  quartz  crystal  clocks,  is  still  much  less 
than  the  cost  of  cesium  clocks.  The  Naval  Observatory  was  asked  to  provide 
information  on  the  predictability  of  such  clocks.  The  initial  experimental 
evaluations  of  such  predictability,  made  by  the  Naval  Observatory  under  OCA 
sponsorship,  indicate  that  useful  predictions  can  be  made  for  a  significant 
period  of  time  [14], 

Timing  stability  and  accuracy  have  great  importance  to  the  survivability 
of  digital  communications.  Communications  system  timing  stability  can  be 
enhanced  by  the  use  of  filtering  if  such  filtering  employs  very  long  time 
constants  (bandwidths  of  a  few  microhertz).  Such  narrow  bandwidths  can  be 
provided  stably  and  economically  in  frequency  control  loops  by  employing 
relatively  low  cost  microprocessors.  For  survivability  of  the  timing 
function,  every  communications  link  between  major  nodes  of  the  network  should 
also  be  capable  of  serving  as  a  timing  reference  link.  To  make  the  most 
effective  use  of  this  capability,  self-organization  should  be  provided. 
Self-organization  will  automatically  adjust  network  parameters  to  compensate 
for  damage  to  portions  of  the  network.  A  most  practical  way  to  provide  this 
control  function  is  by  use  of  a  microprocessor  at  each  major  node.  Since  this 
function  utilizes  an  insignificant  percentage  of  a  microprocessor 's 
capability,  the  microprocessor  can  be  shared  with  other  timing  functions.  It 
can  be  used  for  narrowband  filtering  in  a  phase  control  loop,  and  to  provide  a 
double-ended  feature  which  removes  from  the  timing  coordination  the  time  that 
it  takes  the  signal  to  travel  between  nodes.  This  adds  greatly  to  frequency 
and  phase  stability  of  nodal  clocks  throughout  the  network.  The  double-ended 
feature  also  contributes  greatly  to  the  accuracy  of  the  measurement  of  timing 
errors  at  the  nodes.  The  same  microprocessor  can  be  used  to  make  the  error 
measurements  independent  of  error  corrections  made  at  other  nodes  throughout 
the  network.  This  minimizes  the  propagation  of  timing  errors  through  the 
network.  The  same  microprocessor  can  be  used  for  combining  phase  reference 
(timing)  information  received  over  many  different  paths  in  a  way  which  will 
improve  measurement  accuracy  and  also  simplify  the  reorganization  following 
damage  to  the  network.  Studies  [11]  have  shown  that  all  of  these  features  can 
be  provided  using  only  a  portion  of  a  microprocessor's  capacity.  Therefore, 
the  remaining  capacity  could  be  applied  to  clock  error  prediction  and 
correction.  If  not  enough  remaining  capacity  were  available,  the  cost  of  an 
additional  microprocessor  would  be  very  low  when  compared  with  the  cost  of  a 
cesium  clock.  When  cesium  clocks  are  used,  prediction  techniques  can  also 
improve  their  performance;  but  the  major  cost  saving  advantage  for  the  DCS 
would  be  the  ability  to  use  lower  cost  clocks  in  place  of  cesium  clocks. 

Although  the  Naval  Observatory  has  conducted  a  preliminary  study  under  DCA 
sponsorship  to  determine  clock  predictability,  only  a  limited  number  of  older 
clocks  were  used,  and  there  were  significant  problems  with  some  of  the  clocks 
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during  the  measurements.  No  attempt  was  made  to  optimize  the  parameters  in 
the  ARIMA  model.  Measurements  were  made  once  an  hour  during  that  study  and  no 
attempt  was  made  to  determine  whether  aliasing*  was  having  any  effect  on  the 
result.  It  is  important  to  add  newer  clocks  with  the  latest  technology  to  the 
clocks  which  the  Observatory  already  has.  It  is  also  important  to  take  data 
on  a  larger  number  of  clocks  for  longer  periods  of  time  in  order  to  increase 
confidence  in  the  statistical  studies.  Measurements  should  be  made  on  at 
least  part  of  the  clocks  at  quarter  hour  (or  more  frequent)  intervals  for  use 
in  evaluating  possible  effects  of  aliasing.  Such  an  evaluation  could  be  made 
by  using  every  fourth  (or  less  frequent  measurement)  and  comparing  the 
prediction  results  with  those  obtained  using  every  measurement.  The  effect  of 
averaging  several  samples  prior  to  application  of  the  ARIMA  model  should  also 
be  investigated.  The  work  needs  to  be  further  extended  to  investigate  the 
programming  of  a  microprocessor  to  carry  out  both  clock  error  prediction  and 
the  control  of  a  micro-phase-stepper**  to  correct  the  output  of  the 
free-running  clock  in  a  way  that  could  be  conveniently  applied  to  the  DCS. 


*  Aliasing  —  The  introduction  of  error  in  a  system  employing  discrete 
sampling  of  continuous  data  which  occurs  when  that  data  contains  frequencies 
higher  than  half  the  sampling  rate.  It  occurs  because  the  higher  frequencies 
of  the  original  spectrum,  as  they  occur  around  the  sampling  rate  after 
sampling,  overlap  the  spectrum  of  the  original  continuous  data  and  are 
indistinguishable  from  the  original  frequencies.  Aliasing  can  cause  the 
higher  frequencies  of  the  original  data  to  have  the  same  effect  on  the 
analysis  as  if  they  had  been  low  frequencies  in  the  original  data. 

**Micro-phase-stepper — a  device  which  can  alter  the  phase/time  or  frequency  of 
an  input  reference  signal  by  making  precise  discrete  phase  shifts  which  are 
digitally  controlled.  One  conmercially  available  unit  can  make  shifts  in  a 
5MHz  reference  signal  of  1  picosecond  to  1  microsecond  with  a  resolution  of  1 
picosecond  and  frequency  changes  of  t  1  part  in  10-*7  to  t  1  part  in  10-7 
with  a  resolution  of  1  part  in  lO-*7. 


III.  APPLYING  ARIMA  MODELS  TO  CLOCK  ERROR  PREDICTION 


Autoregressive  Integrated  Moving  Average  (ARIMA)  models  to  be  used  in  this 
program  are  discussed  in  reference  [15],  a  large  text  written  in  the  terms  of 
the  mathematician  and  statistician  and  not  easy  for  most  engineers  to  read.  A 
brief  discussion  of  its  application  to  clock  prediction  written  more  for  the 
engineer  is  provided  in  this  section.  This  explanation  should  be  adequate  to 
give  an  engineer  an  understanding  of  the  application  of  these  models  to  clock 
prediction,  and  it  should  also  put  him  in  a  better  position  to  use  reference 
[15].  Some  of  this  discussion  results  from  the  author's  own  observations  and 
is  not  based  directly  on  the  reference. 

A  model  of  the  form 

Z't  -  tflZ't-1  +  *2z't-2  +  ...  +  dpZ't-p  +  at  where 

Z't  ■  Zfu,  with  Zt  representing  the  value  of  the  process  at  time  t  and 
v  representing  the  mean  about  which  the  process  varies,  is  called  an 
autoregressive  process  of  order  p  because  it  is  a  regression  of  the  variable 
Z'  on  the  p  previous  values  of  itself.  The  0's  are  parameters  of  the  process 
and  at  is  a  random  shock. 

An  autoregressive  operator  of  order  p  may  be  defined  as  0(B)  =  1  -  0iB  - 
02  Bz  -  ..*  -  0pBP  in  which  B  is  a  backward  shift  operator,  i.e., 

BZt  -  Z-t— l »  and  BmZt  =  Zt-m-  Then  the  autoregressive  model  may  be 
written  as  0(B)Z't  -  at 

Similarly,  a  model  of  the  form 

Z't  *  at  -  ©iat-l  -  ©2at-2  -  •••  -  ©qat_q  which  is  dependent 

on  a  finite  number,  q,  of  previous  at's  is  called  a  moving  average  process. 

Using  the  operator  B,  it  may  be  written  as  Z't  =  (l-e(B))at. 

Many  time  series  exhibit  a  nonstationary  behavior  and  do  not  vary  about  a 
fixed  mean.  However,  the  series  formed  by  taking  the  first,  second,  or  higher 
difference  of  these  series  are  frequently  stationary.  Such  behavior  can  be 
represented  by  a  generalized  operator  4(B)  in  which  one  or  more  of  the  zeros 
of  the  polynomial  4(B)  is  unity.  Thus,  4(B)  can  be  written  as  4(B)  = 

0(B) (l-B)d.  Since  the  differences  of  a  discrete  sequence  can  be  related  to 
the  derivative  of  a  corresponding  continuous  function,  the  inverse  of  the 
differences  (sum)  as  sometimes  used  when  working  with  ARIMA  models  can  be 
related  to  an  integral  of  a  corresponding  continuous  function.  In  solving  an 
equation  with  d  differences,  it  is  common  to  insert  initial  values  and  then 
sum  d  times.  Therefore,  a  generalized  autoregressive  function  in  which  d 
zeros  of  the  autoregressive  operator  are  unity  is  said  to  be  integrated  with 
order  d.  A  model  containing  both  a  generalized  autoregressive  part  and  a 
moving  average  part  is  said  to  be  an  Autoregressive  Integrated  Moving  Average 
Model  of  order  p,  d,  q,  i.e.,  ARIMA  (p,  d,  q).  It  has  the  form 
0(B)(l-B)dZ't  *  e(B)at  where  0(B)  is  an  autoregressive  operator  of  order 
p  and  ©(B)  is  a  moving  average  operator  of  order  q. 

Note  that  when  d  >  0  the  equation  contains  only  differences  of  Z't, 
i.e.,  when  d  ■  1  there  are  first  order  differences  corresponding  to 


differences  between  adjacent  values  of  Z't,  and  when  d  >  1  there  are  also 
Signer  order  differences,  such  as  differences  between  first  order  differences, 
etc.  The  significance  of  this  is  that  when  d  >  0,  as  it  always  is  in  clock 
error  measurements,  the  value  of  y  in  Z't  drops  out  and  it  is  only  necessary 
to  be  concerned  with  Zt.  However,  a  value  of  y  can  still  be  subtracted  from 
Zt  if  unreasonably  large  numbers  result  otherwise. 

In  general,  a  moving  average  model  can  be  written  in  terms  of  an 
autoregressive  model  and  vice  versa,  provided  certain  convergence  conditions 
are  met.  To  illustrate  this  take  the  first  order  autoregressive  model 
(l-^lB)Zt  =  at  and  divide  both  sides  by  giving  Zt  = 

( l—^i B )— I at  -  Then  (l-diB)~l  can  be  expanded  into  the  infinite 

sequence  1  +  ^jB+  ^  3g3  +  . ...  Therefore,  the  first 

order  autoregressive  model  can  be  written  as  an  infinite  order  moving  average 

model  with  ©i  =  di,©2  =  i>\  ^,©3  =  <t\  etc.  Similarly,  the 

first  order  moving  average  model  can  be  written  as  an  infinite  order 

autoregressive  model.  This  begins  to  show  the  advantage  of  using  both 

autoregressive  and  moving  average  terms  in  the  model,  an  ARIMA  model.  It  is 

usually  possible  to  get  a  good  model  without  using  a  large  number  of  either 

moving  average  or  autoregressive  terms. 

To  put  things  in  a  form  somewhat  more  familiar  to  electronics  engineers, 
the  ARIMA  model  (l-^(B))Zt  =  (l-©(B))a^  can  be  written  in  the  form  Zt  = 

( (1— ©( B) )/ B) ) )  at  where  the  expression  ( 1  — © ( B ) )/ (l-d( B) )  is  known  as 
the  transfer  function. 

For  continuous  functions  in  electronics  engineering,  the  transfer  function 
is  usually  written  as  H(S)  =  Y(S)/X(S)  where  X(S)  is  the  Laplace  transform  of 
the  input  and  Y(S)  is  the  Laplace  transform  of  the  output.  If  the  poles  of 
H(S),  i.e.,  the  zeros  of  X(S),  have  negative  real  parts,  the  system  is  said  to 
be  stable.  If  their  real  parts  are  positive,  the  system  is  said  to  be 
unstable,  and  the  output  will  increase  without  bound.  With  the  increased  use 
of  digital  technology  in  electronics  engineering,  more  use  is  being  made  of 
discrete  transfer  functions  using  the  Z  transform.  These  are  written  in  the 
form  H(Z)  =  Y(Z)/X(Z)  and  the  transfer  function  of  a  digital  filter  has  the 

M  L 

form  H(Z)  =  (1  ♦  a^ Z- ^  ) / ( 1  +  b^Z^).  Note  that  this  has 

K»|  KM 

the  same  general  form  as  the  transfer  function  of  the  ARIMA  model  if  Z“*  = 

B,  i.e.,  if  the  Z  operator  of  the  Z  transform  and  the  backward  shift  operator, 
B,  are  reciprocals.  Indeed,  this  is  the  case,  because  multiplying  by  Z~1  in 
Z  transforms  corresponds  to  a  unit  step  backward  in  time  just  the  same  as  the 
B  operator  in  the  ARIMA  process.  The  entire  left  side  of  the  s-plane  used  for 
pole  locations  in  continuous  filters  (the  part  where  poles  of  H(S)  can  be 
located  for  a  stable  system)  maps  into  the  inside  of  the  u<iit  circle  in  the 
Z-plane.  For  stability,  the  poles  of  H(Z)  must  fall  inside  the  unit  circle. 
Since  B  is  the  reciprocal  of  Z  (the  Z  transform  operator),  the  zeros  of 
(1-$(B))  must  fall  outside  the  unit  circle  for  stability. 
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From  the  preceding  discussion,  it  can  be  seen  that  the  ARI'MA  model  can  be 
considered  to  represent  the  output,  Z-t,  from  a  linear  filter  whose  input  is 
white  noise,  at-  If  d  =  0  and  all  of  the  zeros  of  (l-tf(B))  are  outside  the 
unit  circle,  the  filter  is  stable;  if  the  zeros  are  inside  the  unit  circle,  it 
is  unstable.  (Alternatively,  it  could  be  considered  to  be  a  process  for 
transforming  Zt  to  white  noise.)  If  d  «  0  and  all  of  the  zeros  of  (l-tf(B)) 
are  outside  the  unit  circle,  the  output,  Zt,  for  a  white  noise  input,  at, 
with  zero  mean  will  be  stationary,  i.e.,  its  statistics  will  be  independent  of 
time.  If  d  =  1,  i.e.,  there  is  one  zero  on  the  unit  circle,  the  first 
difference  of  the  output  Zt  will  be  stationary  with  zero  mean,  but  the  level 
of  Zt  can  be  nonstationary  and  wander  around.  Similarly,  if  d  =  2,  both  the 
level  and  the  slope  are  free  to  wander  around.  The  applicability  of  such  a 
model  to  clock  analysis  and  prediction  begins  to  take  shape. 

The  ARIMA  model  is  a  difference  equation.  Difference  equations  are 
relatives  of  differential  equations  which  are  more  familiar  to  many 
engineers.  Like  differential  equations,  difference  equations  have  a  general 
solution  that  is  the  sum  of  two  parts:  a  complementary  function,  and  a 
particular  integral.  The  complementary  function  is  the  solution  of  the 
equation  when  all  at  are  zero  so  that  only  the  autoregressive  terms  play  a 
part  in  the  complementary  function.  A  distinct  real  root  of  4(B)  *  0  will 
contribute  a  negative  exponential  term  to  the  complementary  function.  A  pair 
of  complex  roots  contributes  a  damped  sine  wave.  Multiple  equal  roots 
contribute  the  product  of  a  polynomial  and  an  exponential.  If  d  equal  roots 
are  unity,  as  d  unity  roots  in  an  integrated  autoregressive  model,  the 
exponent  of  the  exponential  is  zero;  i.e.,  the  exponential  is  a  constant,  1, 
and  the  roots  contribute  a  polynomial  of  degree  d-1.  The  contributions  of 
various  types  of  roots  to  the  complementary  function  are  very  important  to  the 
use  of  the  ARIMA  model  in  predicting  clock  performance.  They  are  important 
because  during  the  period  that  is  being  predicted,  the  inputs,  at,  will  be 
unknown  and  it  will  be  necessary  to  assume  that  they  are  zero.  Therefore,  the 
prediction  will  be  made  using  the  complementary  function  and  a  satisfactory 
choice  of  initial  conditions'.  It  is  not  reasonable  to  expect  any  clock  errors 
that  exist  prior  to  the  start  of  the  prediction  period  to  exponentially 
decrease  during  the  prediction  period.  Similarly,  unless  the  clock  is 
subjected  to  some  form  of  cyclic  environmental  conditions,  errors  would  not  be 
expected  to  have  sine  wave  terms,  either  damped  or  undamped.  Therefore,  it 
would  be  expected  that  clock  errors  could  best  be  approximated  by  some  firm  of 
polynomial,  and  it  would  further  be  expected  that  any  polynomial  that  could  be 
estimated  with  sufficient  accuracy  to  be  useful  would  be  relatively  simple. 

From  the  foregoing  it  can  reasonably  be  expected  that  the  ARIMA  model  will 
have  no  autoregressive  factors  other  than  the  zeros  at  1.  The  ARIMA  model 
should  be  an  ARIMA  (0,d,q),  i.e.,  a  model  with  p  =  0  and  both  d  and  q 
relatively  small  integers.  If  d  =  1,  the  complementary  function  would  be  a 
zero  order  polynomial  which  would  predict  a  constant  phase  error.  If  d  =  2, 
the  complementary  function  would  be  a  first  order  polynomial,  and  it  would 
predict  a  linear  change  in  phase  error  added  to  the  initial  phase  error,  i.e., 
an  initial  phase  error  plus  a  constant  frequency  error.  If  d  =  3,  the 
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complementary  function  would  be  a  second  degree  polynomial,  and  It  would 
predict  a  linear  change  in  frequency  added  to  the  initial  frequency  error  and 
the  initial  phase  error. 

Since  quartz-crystal  clocks  and  rubidium  gas  cells  are  known  to  have 
frequency  drift,  it  would  seem  that  an  ARIMA  (0,3,q)  model  might  be 
appropriate  to  try  for  these  types  of  clocks.  On  the  other  hand,  cesium  beam 
clocks  have  very  little  frequency  drift  so  that  ?  (0,2,q)  model  might  be  more 
appropriate  for  them. 

If  the  nature  of  the  prediction  function  is  determined  entirely  by  the 
complementary  function,  which  in  turn  is  determined  entirely  by  the 
autoregressive  function  of  order  p+d,  where  does  the  moving  average  function 
play  its  part?  The  answer  is,  for  determining  the  initial  values,  i.e.,  how 
the  function  determined  by  the  autoregressive  terms  is  fitted  to  the  measured 
clock  error  data.  For  this  discussion,  let  7^+f  be  the  prediction  for  the 
error  that  will  exist  f  measurement  periods  beyond  the  measurement  time,  t,  at 
which  the  prediction  is  made.  Since  there  have  been  no  measurements  made 
beyond  time  t,  the  values  of  at  must  be  assumed  to  be  zero  after  time  t.  In 
operator  form,  the  ARIMA  (0,3,3)  model  is 
(1  -  Be)3zt  =  (1  -  ©iB  +  ©2B2  -  e3B3)at. 

Written  in  expanded  form,  this  is 

zt  ~  3zt-l  +  3zt-2  “  zt-3  *  at  -  »lat-l  ~  e2at-2  - 
©3at_3 

and  for  prediction  this  can  be  written  as 

Zt+f  -  3Zt+f_i  +  3Zt+f_2  -  Zt+f_3  *  at+f  -  ©iat+f-1  - 
&2at+f-2  “  e3at+l-3 

Then  assuming  that  for  all  values  of  the  subscript  greater  than  t,  a?=0,  and 
zs  “  ^s  (a  predicted  value),  the  following  predicted  values  are  obtained. 

^t+1  *  3zt  -  3zt-l  +  zt-2  "  *lat  ~  ®2at-l  "  »3at-2 

A  A 

zt+2  *  3Zt+i  -  3Zt  +Zt_l  -  ©2at  “  ®3at-l 

AAA 

zt+3  *  3Zf*-2  -  3Zt+l  +  zt  "  ®3at 

A  A  A  A 

Zt+4  *  3Zt+3  -  3Zt+2  +  zt+l 

A  A  A  A 

Zt+f  «  3Zt+f_l  -  3Zt +f-2  +  zt+f-3  for  f>3 

A 

Note  that  all  predictions  beyond  Z^+3  are  formed  completely  from  previous 
predictions  and  no  measured  data  are  directly  included.  Hence,  from  this 
point  forward  the  initial  conditions  have  been  established  and  the  form 
determined  by  the  complementary  function  is  followed;  the  information  that  p  * 
0,  d  *  q  «  3  along  with  the  values  of  ?t+l»  ?t+2»  and  ^t+3  completely 
determine  all  future  predictions.  The  value$Aof  ej,  eo,  and  ©3  are 
important  in  determining  the  values  of  ?t+l>  zt+2  and  ^t+3*  IT  ©3  is 
zero,  the  predicted  errors  for  all  f  >2  are  dependent  on  the  values  of  Z^, 
Tt+i ,  and  Tt+2.  If  both  ©2  and  ©3  are  zero,  the  predicted  errors  for 
all  f  >1  are  dependent  on  the  values  of  Zt_i,  Z^,  and  7t+i.  If  ©j. 


©2,  and  ©3  are  all  zero,  the  predicted  errors  for  ail  f  >  0  are  dependent 
only  on  the  values  of  Zt_2»  Z-t_i *  and  Zt.  It  is  clear  that  if  all 
values  of  ©  are  zero,  any  measurements  made  prior  to  Zt_2  have  no  influence 
on  the  predictions.  Therefore,  it  appears  that  the  values  of  ©  determine  how 
the  measurements  made  prior  to  Zt«2  influence  the  predictions.  In  order  to 
get  an  idea  of  how  ©  can  bring  the  previous  measurements  into  the  prediction 
process,  consider  the  simple  moving  average  model  Zt  *  (l-©iB)at. 

Remember  that  in  making  the  measurements  only  the  actual  clock  errors,  Zt, 
are  observed.  There  is  no  direct  observation  of  the  random  disturbances, 
at-  These  disturbances  must  be  indirectly  inferred  from  the  measurements. 
Solving  the  equation,  Zt  *  H  “  ®lat-l»  simple  first  order 

moving  average  model  to  get  at  gives 
at  *  Zt  +  »iat_i. 

Repeated  substitution  of  this  expression  into  itself  gives: 
at  -  zt  +  ©l(zt-l  +  elat-2 ) 

-  Zt  +  »l ( Zt_l  +  ©1 ( Zt-2  +  ©lat-3)) 

=  Zt  +  ©1 ( Zt— 1  +  ©l(Zt-2  +  ©l(zt-3  +  ©lat-4))) 

*  Zt  +  ©i Zt_l  +  ©1  2Zt_2  +  ©1  3Zt_3  +  “•  +  ©1  nat_n 

In  each  of  these  equations,  all  terms  except  the  last  contain  a  Z,  while  the 
last  term  always  has  an  a  instead  of  Z.  If  n  is  large  and  ©,  is  small,  the 
effect  of  assuming  that  at_n  =  0  can  be  very  small  because  it  is  multiplied 
by  ©i  n.  For  example,  if  ©j  =  0.5  and  n  «  10,  the  error  in  at  from 

assuming  that  at_n  =  0  is  less  than  0.1  percent;  and  if  n  =  20,  it  is  less 

than  0.0001  percent.  However,  the  value  of  at  obtained  in  this  manner 
contains  weighted  values  of  all  Zt  that  occur  after  that  value  of  at_n 
that  was  assumed  to  be  zero.  In  effect,  ©,  determines  the  rate  of  decay  of 
influence  from  older  measurements  in  evaluating  the  latest  at- 

Taking  the  expression  for  the  simple  moving  average  model  ARIMA  (0,0,1) 
and  substituting  for  each  Zt  the  corresponding  expression  for  the 
autoregressive  part  of  an  ARIMA  (0,3,1)  model,  the  value  of  at  can  be 
evaluated  for  that  model  as 

at  =  zt  +  ( ©1 — 3 )*t-i  +  (®1  2  ~  3©1  +  3)zt-2+(®l  3  -  3el 

+3©i  -  l)Zt_3  ♦  ©I f ©1  3  -  3©i  3  +  3©i  -  l)Zt-4  +  •••  + 

©1  "~3(©i  3  -  3©j  c  +  3©i  -  1)  at_n - 

As  with  the  simple  moving  average  model,  ARIMA  (0,0,1),  the  ARIMA  (0,3,1) 
model  has  very  little  error  introduced  in  the  value  of  at  from  assuming 
at-n  *  0,  provided  n  is  large  and  ©j  is  small. 

Reference  [15]  describes  a  method  of  backward  estimation  whereby  measured 
values  are  used  to  estimate  what  the  values  would  have  been  prior  to  the  first 
measurement.  By  using  this  approach,  it  is  possible  to  arrive  at  more 
accurate  values  of  at  more  quickly.  However,  if  a  microprocessor  is 
employed  in  making  the  measurements,  there  should  be  no  problem  in  obtaining 
enough  measurements  for  accurate  estimates  of  at  without  the  use  of  the  much 
more  complicated  procedure  using  backward  estimates. 
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Returning  to  the  discussion  of  the  ARIMA  (0,3,3)  prediction  model,  recall 
that  it  was  observed  that  when  ©2  *  &3  *  0  an  ARIMA  (0,3,1)  model  is 
obtained.  Also,  when  ©2  ■  ©3  ■  0,  all  predictions  after  the  first  are 
based  on  the  last  two  measurements  plus  the  first  prediction,  with  the  first 
prediction  based  on  the  last  three  measurements  and  a  computed  value  of  the 
random  shock,  a*.  This  indicates  the  importance  of  having  accurate 
measurements.  The  first  prediction  gives  the  initial  phase  error.  The  first 
prediction  and  the  last  measurement  give  the  initial  frequency  error.  The 
first  prediction  and  the  last  two  measurements  combine  to  determine  the  rate 
of  frequency  drift.  As  an  example,  assume  that  an  ARIMA  (0,3,1)  prediction 
model  has  an  error  of  size  a  in  the  last  measurement,  Z*.  In  the  evaluation 
of  the  first  prediction  zt+i-in  this  model,  Zt  is  multiplied  by  3 
introducing  a  3a  error  into  Zt+1»  but  Zt  also  appears  in  at  which  is 
multiplied  by  -©2  so  that  the  resulting  error  in  the  prediction  of  the 
initial  phase  error  is  (3-©i)a.  The  initial  prediction  of  the  frequency 
error  comes  from  7t+i  -  Zt.  The  error  in  this  initial  prediction,  which 
results  from  making  an  error,  a,  in  the  last  measurement  Zt,  is  (3  -  ©i)a 
-  A  ■  (2  -  ©i)a.  This  error  in  frequency  prediction  would  cause  an  error  in 
clock  prediction  which  changes  linearly  with  time.  Perhaps  of  greater 
importance,  the  prediction  of  the  frequency  drift  comes  from  the  second 
difference  of  Zt_i,  Zt,  and  ?t+l •  F°r  a  a  error  in  the  last 
measurement,  the  error  in  drift  is  (1-©i)a.  This  error  in  the  prediction  of 
the  drift  would  cause  a  prediction  error  that  would  change  as  the  square  of 
time.  Hence,  it  is  very  important  to  keep  measurement  errors  small.  One 
obvious  method  of  doing  this  in  a  practical  clock  prediction  application  using 
automated  measurements  is  to  make  measurements  once  per  second  and  average 
them  for  several  minutes  before  using  them  as  a  measurement  input  to  the  ARIMA 
model.  This  could  be  particularly  important  for  an  application  such  as  the 
DCS  where  the  measurements  themselves  might  sometimes  be  contaminated  by 
noise.  In  using  this  average,  care  should  be  taken  to  eliminate  from  the 
average  any  individual  measurements  that  are  obviously  displaced  from  their 
neighbors  enough  to  indicate  probable  contamination  of  that  particular 
measurement. 

The  ARIMA  (0,3,3)  model  was  selected  as  the  basis  for  discussion  and  for 
use  in  examples  partly  because  it  seemed  to  have  characteristics  that  would 
make  it  very  useful  for  clock  error  prediction.  It  also  has  the  advantage 
that  by  letting  ©3  -  0  in  the  equations  for  an  ARIMA  (0,3,3)  model,  the 
ARIMA  (0,3,2)  model  is  obtained,  and  by  letting  ©2  =  ©3  =  0,  an  ARIMA 
(0,3,1)  model  is  obtained.  The  ARIMA  (0,0,1)  model  and  the  ARIMA  (0,3,1) 
model  were  chosen  for  the  particular  examples  of  evaluations  of  a*  because 
of  their  simplicity.  They  do  illustrate  that  the  values  of  the  e's  control 
how  earlier  clock  error  measurements  contribute  to  future  predictions  of  clock 
errors.  If  the  values  of  the  random  at's  were  known  and  the  model  were 
absolutely  correct,  the  ARIMA  model  would  be  completely  deterministic  and  the 
prediction  would  be  perfect.  If  the  model  were  absolutely  correct  and  the 
values  of  the  at's  were  accurately  determined  through  the  last  measurement, 
then  the  prediction  would  start  with  the  correct  initial  conditions.  The 
prediction  would  be  correct  except  for  those  random  perturbations  caused  by 
the  at's  that  occur  after  the  beginning  of  the  prediction  period,  and 
therefore  cannot  be  predicted,  the  ARIMA  model  separates  the  deterministic 
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predictable  performance  of  the  clock  from  the  random  unpredictable 
performance.  One  object  of  this  study  program  is  to  determine  the  magnitude 
of  the  unpredictable  part.  Since  this  is  dependent  on  the  at's,  the  rms 
value  of  the  at's  for  a  particular  clock  is  of  interest.  However,  it  is 
also  of  interest  to  make  predictions  of  clock  errors  over  a  long  period  of 
time  and  compare  those  predicted  error  values  with  those  that  do  occur  for  a 
large  number  of  samples. 

The  at's  can  be  obtained  with  a  precision  that  increases  during  a  period 
of  time  during  which  clock  error  measurements  are  made  and  which  will  be 
called  the  calibration  period.  For  the  application  of  clock  prediction 
techniques  to  the  DCS,  the  calibration  period  can  be  the  entire  period  of  time 
during  which  clock  error  measurements  are  made  prior  to  entering  a 
free-running  clock  mode  of  operation. 

The  conditional  expectation  estimate  of  $t+f»  given  all  Z's  up  to  time 
t,  is  the  minimum  mean  square  prediction  at  time  t  for  the  value  at  Zt+f . 

For  a  one  step  ahead  estimate  the  error  is  the  difference  between  the  new 
measured  value  and  its  prediction  made  just  prior  to  the  measurement,  and  this 
error  is  equal  to  the  random  impulse  a++i>  which  is  included  In  the  actual 
measurement  but  which  could  not  be  included  in  the  prediction.  One  method  of 
evaluating  the  at's  is  to  subtract  from  each  new  measurement  the  value  that 
was  predicted  for  it  the  previous  time  period,  at+i  *  Zt+i  -  ?t+l*  This 
process  can  be  started  at  the  initial  measurement,  t  *  0,  by  assuming  all 
unknown  at  information  is  zero.  The  accuracy  of  the  estimates  of  at  will 
increase  as  additional  measurements  are  used.  This  is  illustrated  for  the 
ARIMA  (0,3,3)  model  in  the  following  sequence  of  equations. 

*1+1  -  3Z1 
a2  »  Z2  -  $1+1 
1z*l  .  3Z2  -  3Z]  - 

A 

a3  *  z3  “  z2+l 
z3+l  *  3Z3  “  3Z2  + 

A 

a4  -  z4  ~  z3+l 

^4+1  -  3Z4  -  3Z3  + 

A 

a5  -  z5  *  z4+l 

Z5+1  -  3Z5  -  3Z4  + 

A 

an  ■  Zn  -  Z(n_i)+i 
A 

Zn+1  *  3Zn  •  3zn-l 


©la2 

h  ~  ®la3  "  92a2 

z2  ~  »la4  "  *2a3  -  *3a2 

z3  -  *la5  *  »2a 4  ~  ®3a3 

+  zn-2  -  »lan  ~  ®2an-l  -  »3an-2* 
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As  n  grows  large  the  accuracy  of  the  evaluations  of  at  increases  rapidly, 
provided  the  e's  aren't  too  large. 

The  discussion  to  this  point  should  have  provided  the  reader  with  an 
understanding  of  how  the  ARIMA  model  might  be  used  to  predict  clock  errors, 
but  has  provided  no  indication  of  how  the  model  can  be  optimized  for  a 
particular  clock  or  type  of  clocks.  The  optimization  of  the  ARIMA  model 
comprises  selecting  the  best  set  of  integers  for  the  parameters  p,  d,  and  q, 
and  also  selecting  the  best  values  for  each  and  in  the  model. 

Reasoning  as  to  what  might  be  expected  for  clock  error  performance  has 
indicated  that  for  quartz-crystal  or  rubidium  clocks  a  good  model  might  have 
p  »  0,  d  ■  3,  and  q  ■  some  small  integer.  Similar  reasoning  indicates  that  a 
model  with  p  -  0,  d  *  2,  and  q  *  some  small  integer  might  be  a  good  model  for 
cesium  beam  clocks.  This  needs  to  be  verified  by  determining  how  well  both  of 
these  models  fit  the  measured  data.  In  some  cases,  it  right  be  desirable  to 
compare  results  using  each  of  these  models  with  those  having  small  variations 
in  the  values  of  p,  d,  and  q.  If  the  above  reasoning  about  the  values  of  p 
and  d  are  correct,  only  q  remains  to  be  determined.  Taking  the  third 
differences  of  the  measurements  on  the  clocks  should  provide  a  moving  average 
process.  According  to  reference  [15],  the  autocorrelation  function,  |<,  of 
a  moving  average  process  of  order  q  has  a  finite  value  for  q  <  k  but  is  zero 
for  k  >  q.  Thus,  by  taking  the  autocorrelation  function  of  tFe  third 
differences  (or  second  differences  in  the  case  of  cesium  clocks),  the  value  of 
q  should  be  determined. 

Once  values  of  p,  d,  and  q  have  been  selected,  there  still  remain  some 
unknown  parameters  —  the  values  of  the  tf's,  the  values  of  the  e's  and  ca, 
where  oa  is  the  standard  deviation  of  the  at's.  Reference  [15]  provides  a 
general  discussion  of  methods  for  estimating  these  parameters  using  both 
likelihood  and  Bayesian  approaches.  The  reference  gives  the  unconditional 
log-likelihood  function  as 


l(tf.e,°a)  -  f(d,e)  -  n  lnaa~  S(d,e) 

where  the  i  and  e  in  the  equation  are 
i)  *  >  ^2  * ^3  *  •  •  • » 


®  1 * ®2  *  *  •  •  •  *  &n  ■ 

n 


s<«.»)  .f—  recat  i  *,*.^8] 

2(7  t*-» 


where^d  is  the  set  of  d'th  order  differences  of  observed  data  Z,  and 
EL  At  is  the  expected  value  of  at  conditioned  on  4,  e,  and 

^Z.  The  S(d,e)/2oa  ‘  is  called  the  unconditional  sum  of  squares  and  is 
also  a  function  of  the  d's  and  e's. 


13 


The  model  is  optimized  by  selecting  those  values  of  6  and  ©  which  for  the 
particular  data, maximize  the  log-likelihood  function.  In  computing 
the  unconditional  sum  of  squares,  backward  prediction  is  used  to  determine  the 
values  ofydZ  prior  to  t-0,  i.e.,  t  -  -1,-2, •••.  These  are  used  in  a 
backward  recursion  to  compute  a_j,  a_2,  a_3,***.  These  backward 
values  of  a*  and^Z  can  then  be  used  for  starting  values  in  a  forward 
recursion.  The  unconditional  sum  of  squares  is  obtained  by  summing  the 
squares  of  all  computed  a^'s.  Strictly  speaking,  the  unconditional 
likelihood  is  needed  for  parameter  estimation.  However,  the  term  f($,e)  in 
the  unconditional  log-likelihood  function  is  usually  negligible  for  large 
values  of  n.  Similarly,  if  n  is  large  and  there  are  no  zeros  of  near 
the  unit  circle,  a  conditional  sum  of  squares  can  be  used  which  is  conditional 
upon  the  starting  values  of^Z  and  at.  For  clock  prediction,  both  of 
these  conditions  will  usually  be  satisfied  so  that  the  much  simpler 
conditional  log-likelihood  function  can  be  used.  This  does  not  require  any 
backward  recursive  computations,  and  the  at's  can  be  computed  by  the  method 
already  discussed.  By  summing  the  squares  of  at's  obtained  in  that  way,  the 
conditional  sum  of  squares  is  obtained.  By  making  evaluations  of  the  sum  of 
squares  for  various  values  of  ©j,  ©2,  and  ©3,  the  sum  of  squares 
function  can  be  observed.  The  minimum  sum  of  squares  is  a  close  approximation 
to  the  maximum  likelihood  estimate  of  the  parameters  ©j,  ©2,  and  ©3. 

This  makes  good  intuitive  sense  because  the  at's  represent  the  random 
unpredictable  part  of  the  ARIMA  model.  Surely,  if  another  set  of  parameters 
could  be  selected  which  would  provide  a  lower  computed  value  for  the  sum  of 
the  squares  of  the  at's,  there  would  be  a  better  set  of  parameters  for  the 
model.  Therefore,  the  optimum  set  of  parameters  would  seem  to  be  that  set  of 
parameters  which  minimizes  the  sum  of  the  squares  of  the  computed  at's. 

The  ARIMA  models  offer  the  advantage  for  clock  prediction  that  they  permit 
the  prediction  to  be  accomplished  with  relatively  simple  recursive 
computations. 

As  any  engineer  experienced  with  digital  filtering  techniques  is  well 
aware,  unless  the  sampling  rate  (clock  error  measurement  rate  for  the  clock 
prediction  application),  is  higher  by  at  least  a  factor  of  two  than  the 
highest  frequency  component  of  the  sampled  data,  serious  aliasing  can  result. 
This  aliasing  translates  frequency  components  higher  than  half  the  sampling 
rate  into  lower  frequency  components  where  it  is  very  difficult  and  sometimes 
impossible  to  distinguish  them  from  the  lower  frequency  components  that 
legitimately  belong  there.  Aliasing  is  usually  avoided  either  by  filtering 
away  the  higher  frequencies  prior  to  taking  the  digital  samples  or  selecting 
an  adequately  high  sample  rate.  For  the  application  of  clock  prediction 
techniques  to  the  DCS,  it  is  probably  much  easier  to  use  an  adequately  high 
sample  rate.  Samples  of  once  per  second  or  even  more  frequent  if  necessary 
would  cause  no  particular  problem.  Many  samples  could  be  averaged  together  to 
produce  a  single  sample  for  use  in  the  ARIMA  model,  since  the  averaging  tends 
to  filter  off  higher  frequencies.  Because  there  has  been  no  evaluation  of  the 
part  aliasing  might  play  in  the  ARIMA  model  prediction  of  clock  errors,  the 
possibility  of  its  existence  should  be  kept  in  mind  until  it  is  shown  to  not 
be  a  factor. 


IV.  DISCUSSION  OF  THE  INITIAL  U.S.  NAVAL  OBSERVATORY  STUDY 

Reference  [14]  describes  a  preliminary  evaluation  of  "Predictability  of 
Quartz-Crystal  Oscillators  and  Other  Devices,"  which  was  performed  by  the 
Naval  Observatory  under  sponsorship  of  the  DCA.  This  study  resulted  directly 
from  a  recognition  that  a  communications  timing  system  that  has  the  stability 
and  accuracy  to  adequately  support  the  survivability  requirements  of  the 
Defense  Communications  System  could  also  have  the  capability  to  accurately 
determine  the  errors  in  the  nodal  clocks  during  those  periods  of  time  when  the 
timing  system  is  operating  properly.  It  appeared  that  it  might  be  possible  to 
use  this  information  to  predict  the  errors  that  would  occur  in  those  clocks 
when  it  was  necessary  for  them  to  free-run.  These  predicted  errors  could  then 
be  removed  to  provide  much  better  performance  during  the  free-running  period. 
An  important  question  was  whether  this  increase  In  performance  would  be  great 
enough  to  permit  significant  cost  savings  by  permitting  quartz-crystal  or 
rubidium  clocks  to  be  used  in  place  of  cesium  clocks.  Another  important 
question  was  whether  communications  survivability  would  be  increased  enough  by 
employing  such  clock  error  prediction  techniques  to  justify  their  use  even  if 
they  did  not  permit  use  of  lower  priced  clocks.  The  preliminary  study  seems 
to  indicate  that  either  of  these  reasons  would  probably  justify  use  of  the 
techniques.  However,  much  additional  information  is  still  needed. 

In  the  preliminary  study,  five  different  methods  were  used  to  make  the 
predictions:  (1)  a  1st  degree  polynomial  was  fitted  to  the  data  existing 
prior  to  the  beginning  of  the  free-running  period  and  an  extrapolation  of  this 
fitted  polynomial  was  used  for  the  predicted  values  of  error  during  the 
free-running  period,  (2)  the  same  procedure  was  used  for  a  2nd  degree 
polynomial,  (3)  the  same  procedure  was  used  for  a  3rd  degree  polynomial,  (4) 
the  same  procedure  was  used  for  a  4th  degree  polynomial,  (5)  an  ARIMA  model 
was  used.  An  ARIMA  (0,2,1)  model  was  used  and  no  attempt  was  made  to  optimize 
the  model.  Although  this  is  probably  a  good  model  for  cesium  clocks,  the 
discussion  in  the  preceding  section  indicates  that  an  ARIMA  ( 0,3 ,q )  model 
could  be  a  more  appropriate  choice  for  quartz-crystal  and  rubidium  clocks.  No 
attempt  was  made  to  determine  whether  q  ■  2  or  q  ■  3  might  be  better  than  q  • 
1,  and  no  attempt  was  made  to  optimize  ©i  for  the  various  classes  of 
clocks.  The  value  of  ej  »  0.75  was  used  in  the  ARIMA  model  for  all  clocks. 

For  each  of  the  five  methods  of  making  the  predictions,  a  specific 
prediction  error  was  determined  by  subtracting  the  actual  error,  as  measured 
after  the  prediction  time,  from  the  specific  prediction  obtained  using  the 
model.  Different  sets  of  data,  taken  on  the  same  clock  at  different  times, 
for  the  same  length  of  prediction,  were  used  for  separate  calibration 
periods.  Several  prediction  errors,  determined  for  a  given  clock  with  a  given 
calibration  period  and  a  given  length  of  prediction,  were  used  to  compute  a 
single  root-mean-square  (rms)  error  value  for  each  set  of  parameters. 

The  choice  of  an  ARIMA  (0,3,q)  model  for  either  quartz-crystal  or  rubidium 
clocks  Is  supported  by  this  preliminary  study  even  though  only  the  ARIMA 
(0,2,1)  model  was  employed.  This  support  comes  from  comparing  the  results  of 
the  ARIMA  (0,2,1)  model  with  those  obtained  by  fitting  1st,  2nd,  3rd,  and  4th 
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degree  polynomials  to  the  data  taken  during  calibration  periods.  In  many  of 
the  quartz-crystal  clock  evaluations,  the  2nd  degree  polynomial  prediction 
gave  smaller  errors  than  the  other  polynomials,  and  in  some  cases  smaller  than 
the  ARIMA  (0,2,1)  model.  Since  the  ARIMA  (0,2,1)  model  results  in  only  a  1st 
degree  polynomial,  it  would  seem  that  improvement  could  be  made  if  an  ARIMA 
model  resulting  in  a  2nd  degree  polynomial  were  used.  With  rubidium  clocks, 
the  difference  in  accuracy  of  prediction  between  1st  and  2nd  degree  fitted 
polynomials  was  considerably  less  than  the  corresponding  difference  with 
quartz-crystal  oscillators,  and  the  1st  degree  fitted  polynomial  sometimes 
gave  the  best  prediction.  With  cesium  clocks,  for  longer  prediction  lead 
times,  the  1st  degree  polynomial  frequently  outperformed  the  ARIMA  (0,2,1) 
model.  This  would  seem  to  indicate  either  the  possibility  of  finding  a  better 
ARIMA  model  than  the  ARIMA  (0,2,1)  model  or  the  possibility  of  making  some 
other  improvements  such  as  a  better  value  for  ej.  A  good  ARIMA  model  is 
supposed  to  produce  a  prediction  with  minimum  mean  square  error,  so  it  should 
not  be  possible  to  find  a  model  which  gives  a  lower  mean  square  error  in  its 
prediction. 

This  preliminary  study  of  the  predictabil ity  of  clocks  indicated  clock 
errors  in  high  quality  quartz-crystal  oscillators  can  be  predicted  with 
considerable  accuracy  during  a  free-running  period  following  a  calibration 
period.  There  appears  to  be  a  very  high  probability  that  by  using  these 
predictions  to  remove  errors  during  the  free-running  period,  enough  accuracy 
can  be  obtained  over  a  long  enough  period  of  time  to  permit  the  use  of  these 
clocks  as  an  alternative  to  much  more  expensive  cesium  clocks  for  many 
communications  applications.  Since  the  rubidium  clocks  are  much  more 
predictable  than  high-quality  quart2-crystal  clocks,  there  is  an  even  higher 
probability  that  they  could  be  used  as  an  alternative  to  cesium  clocks  with 
correction  of  the  predicted  errors. 

The  digital  time  error  measurements  for  this  preliminary  evaluation  of  the 
predictability  of  clocks  were  made  once  an  hour.  This  was  assumed  to  be 
frequent  enough  to  produce  satisfactory  results,  but  no  evaluations  were  made 
of  the  possibility  of  significant  aliasing  occurring  of  the  type  discussed  in 
the  preceding  section.  Further  work  is  needed  to  determine  whether  there  is 
significant  aliasing  when  measurements  are  only  made  once  an  hour  during  the 
calibration  period. 

The  preliminary  study  used  a  limited  number  of  older  clocks,  and  there 
were  significant  problems  with  some  of  the  clocks  which  Interrupted 
measurements.  It  is  important  to  add  newer  clocks,  using  the  latest 
technology,  to  the  clocks  which  the  observatory  already  has,  and  to  take  data 
on  a  larger  number  of  clocks  for  longer  periods  of  time  in  order  to  Increase 
confidence  in  the  statistical  studies. 
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V.  THE  PROPOSED  STUDY 


In  addition  to  discussing  the  knowledge  to  be  gained  from  this  study,  this 
section  also  discusses  an  approach  to  getting  the  information.  It  points  out 
some  of  the  things  that  should  be  done  to  acquire  the  desired  information 
economically  in  a  form  that  will  be  useful  to  both  planning  and  design 
engineers. 

The  measurements  needed  to  determine  the  predictability  of  high  quality 
quartz-crystal  oscillators  and  other  devices  should  be  made  relative  to  a 
highly  accurate  reference  such  as  the  Naval  Observatory  master  clock,  which  is 
based  upon  an  ensemble  of  a  large  number  of  highly  accurate  cesium  clocks.  To 
develop  statistical  confidence,  the  measurements  should  be  made  on  a  large 
number  of  clocks.  At  least  part  of  these  clocks  should  employ  the  most  modern 
technology,  and  should  be  of  a  type  that  might  be  employed  in  the  DCS.  The 
measurements  must  be  made  over  a  relatively  long  period  of  time  because  of  the 
statistical  character  of  the  information  to  be  determined.  The  preliminary 
study  by  the  Naval  Observatory  used  clock  error  measurements  made  on  a  limited 
number  of  older  clocks.  A  first  major  step  for  the  proposed  study  program  is 
to  acquire  some  additional  clocks  that  use  the  latest  technology.  Similarly, 
any  additional  test  equipment  needed  to  permit  the  automatic  accumulation  of 
measured  errors  of  each  clock  relative  to  the  Naval  Observatory  master  clock 
should  also  be  acquired.  This  equipment  should  provide  for  making 
measurements  at  frequent  intervals. 

In  order  to  evaluate  possible  adverse  effects  due  to  aliasing  when  data  is 
taken  at  hourly  intervals,  some  data  should  be  taken  at  intervals  of  15 
minutes,  or  more  frequently.  In  order  to  reduce  the  effect  of  errors  in 
making  measurements  of  the  clock  errors,  at  least  some  of  the  measurements 
should  be  made  at  intervals  of  one  minute  or  less.  These  error  measurements 
could  be  averaged  for  a  half  hour  or  an  hour  before  being  applied  to  the  ARIMA 
model.  These  more  frequent  measurements  could  also  be  used  for  a  more 
extensive  evaluation  of  any  adverse  effects  due  to  aliasing.  Another 
advantage  of  these  more  frequent  measurements  is  that  they  also  help  to 
develop  a  procedure  more  directly  related  to  practical  application  in  the 
DCS.  In  the  DCS,  frequent  measurements  are  relatively  convenient  to  make, 
while  the  problem  of  making  errors  in  the  measurement  of  clock  errors  is  much 
more  serious  than  in  a  laboratory  environment.  It  would  be  quite  reasonable 
to  make  clock  error  measurements  once  a  second  and  average  them  for  several 
minutes  before  application  in  an  ARIMA  model.  The  data  needs  to  be  taken  over 
a  long  period  of  time  in  order  to  acquire  enough  data  from  a  limited  number  of 
clocks  to  develop  confidence  in  the  statistical  results. 

In  addition  to  increasing  the  number  of  clocks  to  be  measured  and  in 
addition  to  acquiring  some  test  equipment  to  enhance  the  automatic 
accumulation  of  data,  the  measured  data  needs  to  be  analyzed,  and  plans  for 
the  practical  application  of  this  relatively  new  technology  to  the  DCS  need  to 
be  developed. 
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Analysis  of  the  measured  data  should  concentrate  on  the  application  of 
ARIMA  models.  A  first  step  in  such  an  analysis  should  be  to  develop  a  near 
optimum  ARIMA  model  for  the  data  obtained  from  each  individual  clock.  This 
step  is  needed  so  that  some  evaluation  can  be  made  of  how  much  the  optimum 
model  varies  from  one  clock  to  another.  Do  they  have  the  same  values  for  p, 
d,  and  q?  Are  the  e  parameters  nearly  the  same?  Although  some  comparisons 
should  be  made  between  ARIMA  (0,2,q)  and  ARIMA  (0,3,q)  models,  discussions  in 
a  preceding  section  indicate  that  the  ARIMA  (0,3,q)  model  is  expected  to  be  a 
superior  prediction  model  for  both  quartz-crystal  and  rubidium  clocks,  while 
the  ARIMA  (0,2,q)  is  expected  to  be  superior  for  cesium  clocks.  It  might  also 
be  desirable  to  make  an  evaluation  of  an  ARIMA  (l,2,q)  model  where  one  of  the 
three  zeros  of  the  autoregressive  part  of  the  model  is  not  on  the  unit 
circle.  As  part  of  the  selection  of  the  best  model,  the  value  of  q  needs  to 
be  evaluated.  One  step  that  can  be  taken  in  determining  the  value  of  q  is  to 
take  the  3rd  (or  2nd  if  an  ARIMA  (0,2,q)  model  is  used)  difference  of  the 
measured  clock  error  for  a  particular  clock.  For  this  difference  data, 
compute  the  autocorrelation  function.  For  some  significant  offset  value  of 
the  autocorrelation  function,  its  value,  i.e.,  the  value  of  the 
autocorrelation  function,  should  approach  zero.  The  largest  offset  prior  to 
the  autocorrelation  function  approaching  zero  should  correspond  to  the  value 
of  q  in  the  ARIMA  model. 

Another  method  of  evaluating  the  value  of  q  for  the  model  is  to 
arbitrarily  select  a  number  which  is  equal  to  or  greater  than  the  expected 
optimum  number,  e.g.,  q  *  3,  as  a  starting  point.  Using  this  number  for  q, 
select  the  values  of  the  e's  to  provide  an  optimum  fit  to  the  data.  If  q  is 
chosen  too  large,  some  of  the  values  of  the  e's,  e.g.,  03  and  03,  should 
be  very  close  to  zero.  The  major  problem  with  this  approach  is  that  the 
amount  of  effort  required  to  find  optimum  values  of  the  e's  increases 
tremendously  as  the  number  of  different  e's  increases.  Finding  the  optimum 
value  of  e  when  q  *  1  is  quite  simple  relative  to  finding  the  optimum 
combination  of  values  for  three  different  e's.  Since  the  at's  represent  the 
random  unpredictable  part  of  the  ARIMA  model,  the  e  parameters  which  produce 
the  smallest  root-mean-square  (rms)  value  for  the  at's  must  be  a  very  good 
approximation  to  the  optimum  values.  If  several  values  each  of  ej,  ©2, 
and  63  are  used  in  an  ARIMA  model  with  the  measured  error  data  for  a 
particular  clock,  and  the  rms  values  of  the  a-^'s  are  computed  for  each  set 
of  values,  this  information  can  be  plotted  on  graphs  to  obtain  a  good  estimate 
of  the  optimum  values  of  the  e's.  For  greater  accuracy,  an  additional  set  of 
values  near  the  first  selected  optimum  can  be  used  and  the  process  repeated. 
Reference  [15]  describes  an  iterative  approach  for  computing  optimum  values  of 
the  e's,  and  this  approach  should  be  investigated  for  its  possible  application 
to  this  study. 

Once  an  optimum  ARIMA  model  has  been  selected  for  each  individual  clock, 
these  optimum  models  should  be  compared  to  determine  whether  a  single 
compromise  ARIMA  model  for  all  clocks  of  a  given  type  could  be  expected  to 
provide  near  optimum  results,  i.e.,  one  ARIMA  model  for  all  quartz-crystal 
clocks,  another  for  rubidium  clocks  and  another  for  all  cesium  clocks. 
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The  prediction  errors  for  all  optimum  models  selected  for  individual 
clocks  should  be  determined  by  comparing  predicted  values  with  actual  measured 
data.  To  do  this,  a  sequence  of  measured  clock  error  data  is  used  to  perform 
the  calibration  process,  i.e.,  evaluate  the  at's  for  use  in  determining 
initial  conditions  for  the  predictions.  Following  the  last  measured  clock 
error  datum  used  for  calibration,  prediction  is  started  by  assuming  all  future 
a^'s  to  be  zero  in  the  recursive  process.  The  datum  actually  measured  for 
each  of  the  predicted  values  is  subtracted  from  the  predicted  value  to 
determine  the  error  in  the  prediction.  For  each  clock,  this  should  be  done 
for  several  sets  of  data  taken  at  different  times.  It  should  also  be  done  for 
a  family  of  different  lengths  of  calibration  period.  For  at  least  a  few 
clocks,  the  error  in  prediction  should  be  plotted  as  a  function  of  prediction 
time  in  order  to  provide  additional  insight  as  to  the  character  of  these 
errors,  the  unpredictable  part  of  the  clock  errors.  In  addition  to  plotting 
individual  error  curves,  many  different  sets  of  prediction  errors  for  a 
particular  clock  should  be  combined  into  an  rms  evaluation  of  the  prediction 
error  for  that  clock  as  a  function  of  prediction  time. 

It  is  desirable  to  perform  the  above  evaluations  using  the  optimum  ARIMA 
model  for  each  individual  clock,  and  also  using  the  compromise  ARIMA  model  for 
each  type  of  clock  so  that  the  amount  of  degradation  resulting  from  the 
compromise  can  be  evaluated. 

When  all  of  the  above  measurements  and  computations  have  been  completed, 
enough  information  should  be  available  to  make  estimates  of  the  accuracy  with 
which  errors  of  a  particular  type  of  clock  can  be  predicted.  Knowing,  with 
some  degree  of  confidence,  the  accuracy  with  which  clock  errors  for  a 
particular  type  of  clock  can  be  predicted  should  be  a  very  valuable  tool  for 
design  engineers  considering  the  use  of  such  predictions  in  a  practical 
application.  The  information  is  needed  to  determine  whether  the  predictions 
can  be  made  with  sufficient  accuracy  over  a  long  enough  period  of  time  to 
satisfy  the  needs  of  a  particular  application,  and  to  evaluate  any  degradation 
that  might  result  from  the  substitution  of  lower  cost  clocks  for  expensive 
ones. 

Since  no  evaluations  have  been  made  to  determine  the  effects  of  aliasing 
resulting  from  measuring  clock  errors  once  an  hour,  this  new  evaluation  of 
clock  error  prediction  should  include  determining  whether  such  effects  exist 
and  if  they  do,  their  nature.  One  method  of  doing  this  is  to  make  several 
sequences  of  error  predictions  using  measured  error  data  taken  much  more 
frequently  than  once  an  hour,  then  repeating  each  of  these  predictions  using 
only  a  fraction  of  the  measured  data,  e.g.,  every  5^  or  every  l0th 
measurement.  By  comparing  the  results  of  these  different  sets  of  predictions, 
it  should  be  possible  to  determine  whether  any  degradation  In  the  accuracy  of 
prediction  has  occurred  as  a  result  of  using  the  less  frequent  measurement 
data.  By  using  several  different  rates  of  clock  error  measurements,  i.e.,  by 
allowing  different  numbers  of  unused  measured  data  points  between  those  used 
for  making  predictions,  it  should  be  possible  to  evaluate  any  degradation  due 
to  aliasing  as  a  function  of  the  frequency  with  which  clock  error  measurements 
are  made. 


Evaluations  are  also  needed  to  determine  the  effect  of  perturbations  in 
the  measurement  of  clock  error  data.  The  perturbations  might  arise  from 
errors  produced  by  the  measurement  equipment,  or  as  the  result  of  the  way  the 
measurement  equipment  is  used.  In  a  communications  network,  some  perturbation 
can  be  expected  from  the  noise  on  the  communications  link  used  for  clock 
comparisons.  The  perturbations  due  to  transmission  link  noise  can  be  expected 
to  vary  with  the  transmission  medium  employed.  They  are  expected  to  be  very, 
very  small  on  fiber  optics  cable  transmission  links,  but  could  develop  to  a 
significant  value  on  over-the-horizon  tropospheric  scatter  microwave  links. 
Beyond  determining  how  such  perturbations  affect  the  accuracy  of  the 
prediction  results,  it  is  desirable  to  evaluate  methods  of  minimizing  such 
degradation.  As  mentioned  before,  one  method  is  to  make  very  frequent 
measurements  and  to  average  a  large  number  of  them  together  to  obtain  each 
measurement  data  point  used  in  the  ARIMA  model.  In  addition  to  such 
averaging,  steps  should  be  taken  to  increase  the  robustness  of  the  prediction 
process  by  eliminating  from  the  averages  those  measurements  which  are  greatly 
displaced  from  the  mean  of  the  measurements,  or  by  using  some  other  method  of 
robust  statistical  inference. 

For  applications  of  clock  error  prediction  techniques  in  communications 
systems  such  as  the  DCS,  it  can  be  expected  that  it  will  usually  be  possible 
to  have  an  extended  calibration  period  prior  to  the  free-running  clock  period 
during  which  the  predictions  would  be  used.  This  calibration  period  could 
usually  be  expected  to  last  several  months  or  even  years.  It  might  be 
possible  to  develop  a  relatively  simple  microprocessor  algorithm  to  update  and 
more  precisely  determine  the  optimum  values  of  the  e's  in  the  ARIMA  model  for 
the  particular  clock  on  which  the  predictions  would  be  used.  This  might  make 
it  possible  to  obtain  the  most  optimum  model  practicable  for  use  at  the 
specific  time  when  the  prediction  is  needed.  This  possibility  should  be 
investigated  as  a  part  of  this  proqram. 

In  addition  to  the  evaluations  of  the  predictability  of  the  clocks  as 
discussed  above,  investigations  need  to  be  made  into  the  practical 
applications  of  these  techniques  in  nodes  of  military  communications  networks, 
such  as  the  DCS,  to  aid  in  enhancing  the  survivability  of  the  timing  function 
when  the  network  must  withstand  massive  destruction  due  to  enemy  action. 
Relatively  simple  algorithms  that  can  be  applied  to  simple  reliable  low  cost 
microprocessors  are  needed. 

Investigations  should  also  be  made  into  sharing  the  microprocessor  used 
for  applying  the  ARIMA  prediction  model  to  other  purposes.  Those  other 
purposes  might  include  the  ability  to  assure  stable  and  accurate  timing  from  a 
node's  local  clock  under  all  circumstances  including  a  destructive  attack  upon 
the  network.  In  order  to  make  good  predictions  of  clock  errors  during  a 
free-running  clock  mode  of  operation,  good  error  measurements  must  be 
available  during  a  calibration  period  preceding  the  free-running  period.  This 
implies  that  measurement  errors  due  to  perturbations  in  the  communications 
transmission  medium  should  be  kept  small.  To  accomplish  this,  the  local  clock 
at  each  node  controls  the  timing  of  all  transmitted  synchronization  codes. 


20 


The  arrival  times  of  all  received  synchronization  codes  are  measured  relative 
to  the  local  clock.  The  measurement  includes  the  signal  transit  time  and  any 
difference  between  the  two  clocks.  The  measurement  is  communicated  to  the 
other  end  of  the  link.  Subtracting  the  measurement  made  at  one  end  o*  the 
link  from  that  made  at  the  other  end  and  dividing  by  two  gives  the  measured 
difference  between  the  two  clocks.  If  all  nodes  transmit  their  measured  but 
uncorrected  (known)  errors  to  their  neighbors,  any  node  wishing  to  use  a 
particular  neighbor  for  time  reference  information  can  use  this  information 
from  that  neighbor,  along  with  the  measured  difference  between  its  own  clock 
and  that  of  the  neighbor,  to  determine  its  own  measured  but  uncorrected 
error.  The  microprocessor  should  be  programmed  to  take  part  in  this 
measurement  procedure  to  assure  that  good  measurement  data  are  available  for 
application  in  the  ARIMA  model. 

During  periods  of  normal  operation  when  measured  errors  in  the  local  clock 
can  be  used  to  control  the  output  phase  of  the  local  clock,  i.e.,  when  it  is 
unnecessary  for  the  local  clock  to  free-run,  the  clock  error  measurements  need 
to  be  filtered  before  being  used  for  clock  error  corrections.  The 
microprocessor  should  be  programmed  to  provide  this  filtering  to  assure 
stability  unperturbed  by  the  transmission  media  or  changes  in  other  clocks. 

The  most  stable  method  of  applying  clock  corrections,  i.e.,  the  method 
which  least  disturbs  the  basic  stability  of  the  clock,  is  to  provide  a  phase 
correction  in  tandem  with  the  output  of  the  basic  clock.  Digital  control  of 
micro-phase-steppers  is  a  most  accurate  and  effective  method  of  doing  this. 

The  local  clock  signal  is  normally  taken  to  be  the  output  of  this  micro-phase- 
stepper,  while  the  stability  of  the  local  clock  output  is  determined  by  the 
basic  clock  or  oscillator  preceding  the  micro-phase-stepper.  The  measurement 
of  the  error  in  the  output  of  the  mi cro-phase-stepper  and  the  known  amount  of 
correction  provided  by  it  should  be  used  to  compute  the  error  in  the 
uncorrected  local  clock.  This  evaluation  of  the  uncorrected  local  clock  error 
is  used  during  the  normal  periods  of  operation  for  calibration  of  the  ARIMA 
model  to  be  used  for  making  predictions  of  errors  during  a  free-running  period 
to  follow.  It  is  also  used  for  determining  clock  error  corrections  to  be  made 
by  the  micro-phase-stepper  during  normal  operations.  The  microprocessor 
should  be  programmed  to  perform  these  functions. 

During  the  normal  mode  of  operation,  phase  error  measurement  information 
will  be  received  from  several  different  neighboring  nodes.  Therefore,  the 
accuracy  and  stability  of  the  measurement  information  can  be  enhanced  by 
optimally  combining  this  information  received  over  several  different  paths.  A 
method  for  doing  this  is  discussed  in  references  [4],  [7],  and  [13].  The 
microprocessor  should  be  programmed  to  perform  the  computations  which  provide 
this  enhancement  of  measurement  precision.  When  this  method  of  combining  the 
clock  error  measurement  information  from  different  communications  paths  is 
used,  it  also  provides  possibilities  for  quantitative  timing  system  self¬ 
monitoring  which  can  allow  an  alarm  to  be  given  far  in  advance  of  any  actual 
degradation  of  the  communications  provided  by  the  system.  The  microprocessor 
should  also  be  programmed  to  provide  this  function. 
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In  a  military  communications  system  subject  to  massive  destruction, 
provision  must  be  made  to  automatically  reorganize  the  timing  function  so  that 
a  new  network  master  is  selected  when  necessary,  and  the  transfer  of  timing 
information  through  the  network  is  not  degraded  more  than  necessary.  Methods 
for  doing  this  are  discussed  in  references  [4],  [7],  and  [13].  The 
microprocessor  should  provide  this  function  also. 

The  microprocessor  should  also  be  programmed  to  minimize  the  effect  of 
erroneous  or  perturbed  timing  error  measurements  made  during  normal  periods  of 
operation.  Also  during  normal  periods  of  operation,  it  should  filter  the 
measurement  data  prior  to  its  use  for  controlling  a  micro-phase-stepper  to 
correct  clock  errors.  It  should  control  the  amount  of  correction  applied  by 
the  micro-phase-stepper  during  the  controlled  mode  of  operation.  It  should 
also  calibrate  the  clock  prediction  model  during  the  controlled  mode  of 
operation.  It  should  determine  when  it  is  necessary  to  change  from  the 
controlled  mode  of  operat’on  to  the  free-running  mode  of  operation,  and 
provide  a  smooth  transition  from  the  controlled  mode  to  the  free-running 
mode.  During  the  free-running  mode,  it  should  predict  the  clock  errors  and 
control  the  micro-phase-stepper  to  remove  the  predicted  errors.  It  should 
also  determine  when  to  return  to  the  controlled  mode,  and  should  provide  a 
smooth  transition  to  the  controlled  mode.  It  should  also  provide  an  alarm 
capability  to  attract  attention  to  any  timing  system  faults  which  occur,  and 
provide  automat)''  diagnosis  of  the  faults  for  maintenance  personnel. 

Once  a  capability  to  use  one  or  more  microprocessors,  along  with  required 
supporting  equipment  (such  as  time  interval  counters  and  micro-phase-steppers) 
to  perform  all  of  the  desired  timing  system  functions  has  been  demonstrated,  a 
comprehensive  report  will  be  needed.  This  report  should  present  all  of  the 
information  related  to  measuring  clock  errors  during  normal  operation, 
predicting  clock  errors  and  removing  them  during  free-running  operation,  and 
otherwise  enhancing  the  performance  of  the  timing  system  while  improving  its 
survivability  and  lowering  its  costs.  To  the  extent  possible,  the  report 
should  be  written  in  engineering  terms,  providing  charts  and  tables  with  basic 
parameters  and  relationships;  and  it  should  provide  adequate  guidance  on  their 
use  to  make  it  a  good  tool  for  both  the  planning  and  system  design  engineers. 
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VI.  CONCLUSIONS  AND  RECOMMENDATIONS 

The  survival  of  the  communications  function  of  a  military  synchronous 
switched  digital  network  subjected  to  massive  enemy  destruction  of  parts  of 
the  network  is  dependent  on  the  survivability  of  the  timing  function.  The 
timing  function  survivability  can  be  greatly  enhanced  by  providing  a  backup 
free-running  clock  mode  of  operation  at  each  node  of  the  network.  For  higher 
data  rates  such  as  those  expected  between  major  nodes  of  the  future  DCS,  the 
usefulness  of  a  clock  in  the  free-running  mode  is  highly  dependent  on  its 
stability  and  accuracy  during  that  free-running  period.  This  clock  accuracy 
also  aids  the  rapid  resynchronization  of  spread  spectrum  equipment, 
cryptographic  equipment,  and  other  synchronous  devices  when  a  node  reenters 
the  network  after  temporarily  being  isolated  from  the  network.  Techniques  are 
available  for  accurately  measuring  errors  in  local  clocks  during  a  normal 
controlled  mode  of  timing  system  operation.  These  measurements  can  be  used  to 
calibrate  a  mathematical  model  of  the  error  performance  of  the  free-running 
clock.  During  any  contingency  requiring  a  free-running  clock  following  enemy 
destruction  of  parts  of  the  network  or  following  random  equipment  failures, 
these  predicted  errors  can  be  removed,  thereby  greatly  enhancing  the  timing 
accuracy. 

There  is  a  very  high  probability  that  high  quality  quartz-crystal  or 
rubidium  clocks  could  be  used  as  a  much  lower  cost  alternative  to  expensive 
cesium  beam  clocks  in  many  DCS  applications  when  clock  error  prediction  and 
correction  techniques  are  employed.  This  could  result  in  lower  costs,  and 
enhancement  of  the  short  term  accuracy.  In  those  cases  where  quartz-crystal 
or  rubidium  clocks  cannot  be  used  satisfactorily  because  their  errors  cannot 
be  predicted  with  sufficient  accuracy  for  an  adequately  long  period  of  time, 
the  same  techniques  can  be  used  to  enhance  the  accuracy  of  cesium  clocks. 

An  Autoregressive  Integrated  Moving  Average  (ARIMA)  model  is  a  good  tool 
to  use  for  these  predictions.  When  properly  applied,  it  should  produce  a 
minimum  mean  square  error  prediction  of  the  clock  errors.  It  has  the 
advantage  of  permitting  recursive  calculation  to  be  used  during  both  the 
calibration  period  and  the  prediction  period,  and  there  is  some  possibility  of 
developing  automatic  methods  for  model  parameter  optimization  during  lengthy 
cal ibration  periods. 

Although  these  prediction  techniques  are  extremely  promising  for  a  timing 
system  for  a  digital  DCS,  considerable  information  is  required  before  they  can 
be  applied  most  effectively.  It  is  recommended  that  the  general  program  for 
acquiring  the  information  described  in  this  document  be  initiated  and  carried 
to  a  meaningful  completion. 
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